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Abstract. In the steady Couette flow of a granular gas the sign of the heat flux 
gradient is governed by the competition between viscous heating and inelastic 
cooling. We show from the Boltzmann equation for inelastic Maxwell particles 
that a special class of states exists where the viscous heating and the inelastic 
cooling exactly compensate each other at every point, resulting in a uniform heat 
flux. In this state the (reduced) shear rate is enslaved to the coefficient of resti- 
tution a, so that the only free parameter is the (reduced) thermal gradient e. It 
turns out that the reduced moments of order k are polynomials of degree k — 2 in 
e, with coefficients that are nonlinear functions of a. In particular, the rheological 
properties (k — 2) are independent of e and coincide exactly with those of the 
simple shear flow. The heat flux (k = 3) is linear in the thermal gradient (gen- 
eralized Fourier's law), but with an effective thermal conductivity differing from 
the Navier-Stokes one. In addition, a heat flux component parallel to the flow 
velocity and normal to the thermal gradient exists. The theoretical predictions 
are validated by comparison with direct Monte Carlo simulations for the same 
model. 

1 Introduction 

The influence of collisional dissipation on the dynamical properties of granular matter can be 
modeled by a fluid of inelastic hard spheres (IHS) with a constant coefficient of normal restitu- 
tion a < 1. For sufficiently low densities, the Boltzmann equation (BE) has been conveniently 
generalized PQ to account for the inelasticity of binary collisions and the Navier-Stokes (NS) 
transport coefficients have been obtained in terms of the coefficient of restitution |2I3I4] by 
means of the Chapman-Enskog method. However, while the NS equations (constitutive equa- 
tions that are linear in the hydrodynamic gradients) have been shown to be quite useful for 
describing several problems, in most situations of practical interest in granular gases (such as 
steady states) large gradients occur and more complex constitutive equations are required. This 
does not signal a breakdown of hydrodynamics, only a failure of the NS approximation. 

One of the well-known examples of steady states is the simple or uniform shear flow (USF) 
problem [516] . characterized by a linear velocity field (i.e., du x /dy = const) and constant density 
n and temperature T. This type of steady shear flow with zero heat flux can only occur when 
collisional cooling (which is fixed by the mechanical properties of the particles making up the 
granular gas) is exactly balanced by viscous heating (which is controlled by the shearing). 
Consequently, due to this intrinsic relationship between the shear field and dissipation, the 
(dimcnsionless) strength of the velocity gradient is set by the collisional cooling in the steady 
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Fig. 1. The planar Couette flow is driven by two horizontal plates, separated by a distance h. Both 
act like sources of temperature and shear on a low density granular gas filling the space between them. 



state. This means that the corresponding hydrodynamic steady shear flow state is inherently 
non-Newtonian [5] (and so, beyond the scope of the NS equations) in inelastic granular gases. 

The steady planar Couette flow (cf. Fig. [T|) is much more complex than the USF. In contrast 
to the latter, the temperature is not uniform and thus a heat flux vector q coexists with the 
pressure tensor Pij [?]• In general, inelastic cooling and viscous heating are unbalanced, their 
difference dictating the sign of the divergence of the heat flux |7l8j . More explicitly, the energy 
balance equation reads 



(nT + P xy —^, (1) 



where £ is the inelastic cooling rate and d is the dimensionality of the system. Moreover, 
conservation of momentum implies P xy — const and P yy = const. Note that, while the inelastic 
cooling term (YtT is inherently positive, the viscous heating term P xy du x /dy is inherently 
negative since the sign of the shear stress P xy is opposite to that of the shear rate du x /dy. 

It can be shown that, for appropriate values of the boundary conditions for the velocity, 
both terms on the right-hand side of Eq. (flj can exactly balance at all points of the system 
[5] . This yields a whole new set of steady states with a uniform heat flux and other important 
common properties [9]. For these flows, the streamwise heat flux component q x (absent at NS 
order) is also uniform and in general different from zero. Of course, in the absence of thermal 
gradient, q y = and one recovers the steady state condition defining the USF problem. On 
the other hand, in the elastic limit (a = 1 and £ = 0) but nonzero thermal gradient, the 
conventional Fourier flow description (where du x /dy — 0, P xy = 0, and q x = 0) for an ordinary 
gas is also recovered. Thus, the above two well-known steady states (inelastic USF and elastic 
Fourier flow) can be seen as particular situations of this class of steady Couette flows. Apart 
from the condition q = const, a unifying feature of this new class is that, while neither u x (y) 
nor T(y) are linear, a parametric plot of T vs u x shows a linear relationship [8 9 . For this 
reason, we referred to this class of flows as LTu. The slope of the linear plot T{u x ) goes from 
zero in the inelastic USF limit (constant temperature) to infinity in the elastic Fourier flow 
(zero macroscopic velocity). 

Unfortunately, the analysis of the steady Couette flow from the BE for IHS is quite intricate 
and so evidence on the existence of the LTu class was provided in Ref. [§] by molecular dynam- 
ics simulations, a numerical solution of the BE by means of the direct simulation Monte Carlo 
(DSMC) method [TU], and an approximate solution obtained by Grad's method. Of course, it 
is desirable to support the LTu class by means of analytical results cleanly derived from the 
BE. Like in the elastic case, a way to overcome such difficulties is by considering the so-called 
Maxwell models, i.e., models for which the collision rate is independent of the relative velocity 
of the two colliding particles [1 Ill2ll3ll4ll5ll6ll7ll8ll9l20l21l22l23l24l25j . Thanks to this sim- 
plification, nonlinear transport properties can be exactly obtained in some particular problems 
for elastic |26l27j and inelastic [19121122] gases without the need of introducing additional, and 
sometimes uncontrolled, approximations. Apart from their academic interest, it has been also 
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shown that in some cases the results derived from inelastic Maxwell models (IMM) compare 
well with those obtained from IHS 19 22]. 

In this paper we revisit the steady Couette flow with uniform heat flux (LTu class) in the 
framework of the BE for IMM. Given that the shear rate and the coefficient of restitution 
are related through the condition dq y /dy = in ([1}, the only nonequilibrium free parameter 
in the problem turns out to be the (reduced) thermal gradient e cx dT/dy. Based on previous 
simulation and approximate results obtained for IHS [9 , it is shown that the steady BE for IMM 
admits an exact solution characterized by constant pressure, du x /ds = const, and dT/ds = 
const, where s is a conveniently scaled space variable. An analysis of the infinite hierarchy of 
velocity moments shows that the reduced moments of order k = 2r + £ are polynomials of parity 

1 and degree k — 2 in e, whose coefficients are nonlinear functions of a. The evaluation of the 
coefficients associated with the pressure tensor (2r + £ = 2) and with the heat flux (2r + 1 = 3) 
is one of the main goals of this paper. The knowledge of the second-, third-, and fourth-order 
collisional moments for <i-dimensional IMM [53] allows us to explicitly get those coefficients. In 
particular, the expression for the pressure tensor coincides with the one previously obtained for 
the USF |21j , while the heat flux is linear in the thermal gradient with two nonzero components 
(q y and q x ). The constitutive relation between q y and d y T can be seen as a generalized Fourier's 
law and defines an effective thermal conductivity different from the NS coefficient. The existence 
of a component of the heat flux normal to the thermal gradient is a direct consequence of the 
inherent non-Newtonian features of the Couette flow since this component is absent in the NS 
description [2]. 

The above exact solution of the moment hierarchy resulting from the BE for IMM is free 
from boundary layers and thus it holds in the bulk region of the granular fluid. In order to check 
its reliability under realistic boundary conditions, the BE for IMM has also been numerically 
solved for a three-dimensional system by means of the DSMC method. The comparison between 
theory and simulation allows us to confirm that the analytical results are not mathematical 
artifacts but actually describe the steady state achieved by the system in the bulk domain. As 
will be seen, the good agreement found between both approaches clearly supports the validity 
of the assumptions made when working out the theoretical solution. 

The plan of the paper is as follows. In Section [21 the BE for IMM is introduced and the 
explicit expressions for the collisional moments through third-order are given. Section [3J deals 
with the problem we are interested in, namely the steady Couette flow with uniform heat flux 
(LTu class). It is consistently proven that the hierarchy of moment equations associated with 
the problem admits a solution where the moments are polynomials in the thermal gradient 
with coefficients that are nonlinear functions of the coefficient of restitution. Next, the explicit 
expressions for the momentum and heat fluxes are obtained in Sectional while the comparison 
with Monte Carlo simulations is presented in Section [5] The paper is closed in Section [5] with 
a brief discussion of the results derived here. 

2 The inelastic Maxwell model 

Let us consider a granular fluid modeled as a gas of inelastic Maxwell particles. A constant 
parameter, the coefficient of normal restitution a, accounts for the inelasticity in the collisions. 
Its values range from a = (purely inelastic collision) to a — 1 (purely elastic collision). At a 
kinetic-theory level, all the relevant information on the state of the system is conveyed by the 
one-particle velocity distribution function /(r, v; t). In the low-density regime, the inelastic BE 
for IMM reads 



(9 t +vV)/(r,v;i) = J[v|/,/] ) 



(2) 



where the Boltzmann collision operator J[v|/, /] is given by [11113115116117115] 




(3) 



Here, 




(4) 
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is the number density, v is the collision frequency (assumed to be independent of a), fid = 
27r d / 2 / T(d/2) is the total solid angle in d dimensions, and the primes on the velocities denote 
the initial values {v' 1; v 2 } that lead to {vi,V2} following a binary collision: 

v[ = Vl - - (14- a" 1 ) (ff-g)ff, v 2 = v 2 4- - (l + a" 1 ) (£ • g)& , (5) 

where g = vi — v 2 is the relative velocity of the colliding pair and <r is a unit vector directed 
along the centers of the two colliding particles. 

Apart from n, the basic moments of / are the flow velocity 

u = - I dvv/(v) (6) 



and the granular temperature 

T=£ jdv^v), (7) 

where V = v — u(r) is the peculiar velocity. The momentum and energy fluxes are characterized 
by the pressure tensor 

Pij = mj torV&jf{y) (8) 

and the heat flux 

q=^|dvV 2 V/(v), (9) 
respectively. Finally, the rate of energy dissipated due to collisions defines the cooling rate ( as 



t=-£f j dv7 2 J[v|/,/]. (10) 

The collision frequency v is proportional to density and, in general, is also a function of the 
granular temperature T. In fact, in order to make contact with the IHS model, one must take 
v oc nT 1 / 2 [TT]. In particular, the choice [3U] 

d + 2 4l2 d d _ x [f 

»=—ms, "Hs= ^ (d + 2) ™ \J~, (ii) 

where a is the diameter of the spheres, yields the same expression for the cooling rate as the 
one found for IHS (evaluated in the local equilibrium approximation). This will be the choice 
made here, although most of the results are actually independent of the explicit form of v. The 
collision frequency z^hs is the one associated with the NS shear viscosity of an ordinary gas 
(a = 1) of hard spheres, i.e., ?7ns =p/' y HS- 

The main advantage of the BE for Maxwell models (both elastic and inelastic) is that the 
(collisional) moments of the operator J[f, f] can be exactly evaluated in terms of the moments 
of /, without the explicit knowledge of the latter [35]. In particular, the second- and third-order 
collisional moments are [20123] 



ij dvViVjj\y\f, f] = -v 0l 2(Pij -pSij) - cpSn, 



(12) 



y J dV ViVjVkJiVlf, f] = -^Qijk - V2ll d + l° 13 fe^fc + Ufa + QkSii) , (13) 
where p = nT — ^tr P is the hydrostatic pressure, 

777 f 

Qijk = j J dVViVjVkfiV) (14) 
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is a third-rank tensor, and 

A (l + a) 2 3 A (l + a) 2 (d-l) 3 

" l 2 = C+ 2(^f' ^ 1 = 2 C+ 2d(d + 2) ^I3=r°l 2 - (16) 

In Eq. (|16[) the collision frequencies i^ |2 an d ^211 are decomposed into a part inherent to the 
collisional cooling plus the genuine part associated with the momentum and energy collisional 
transfers. For later use, it is convenient to define 

(1 + q) 2 

Id = — C = — -V 

0,2 ^ 2(d + 2) 
2 

^HS- (17) 

Thus the collision frequency i/rs defined by the first equality of Eq. (1111) is actually the value 
of uj in the elastic limit a — > 1. 

The fourth-order collisional moments were also explicitly evaluated in Ref. |23j but will not 
be displayed here [25] . 



3 Couette flow with uniform heat flux 

The planar Couette flow considered in this paper corresponds to a granular gas enclosed between 
two parallel, infinite plates (normal to the y axis) at y — ±h/2 in relative motion along the x 
direction, and kept at different temperatures (cf. Fig. [I}. The resulting flow velocity is along 
the x axis and, from symmetry, it is expected that the hydrodynamic fields only vary in the 
y direction. Consequently, the velocity distribution function is also assumed to depend on the 
coordinate y only. Moreover, we focus on the steady state, so Eq. @ becomes 



v y d a f=-J[f ) f] 

UJ 



1 2(rf + 2) 



J dv 2 J da [a-'/Kl/K) - /(vx)/(v 2 )] , (18) 



nf2 d (1 + a) 2 

where we have introduced the scaled variable s as 

ds = cudy. (19) 

Note that Eq. ([TBI is "universal" in the sense that it is independent of the precise choice of v. 
Now we assume that an exact solution of Eq. (fTB)) exists characterized by 

du x ~ dT 
p = const, — — = a = const, — — = const. (20) 
os os 

The constant a represents a Knudsen number associated with the shear rate. There is another 
Knudsen number associated with the thermal gradient, namely 

-dlnT 



e=y/2T/m--—. (21) 
os 

This quantity is not constant since d s T = const implies e oc T -1 / 2 . From Eqs. (|2T)|) and (f2"T|) we 
get 

dT eJmT/2 

— — = v - 1 = const. (22) 
ou x a 
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This means that when the spatial variable (y or s) is eliminated to express T as a function of 
u x one gets a linear relationship. Therefore, the solution to the BE (Tl8|) consistent with the 
hydrodynamic profiles (|2"0)) defines the LTu class of Couette flows described in the Introduction. 

As we will see, the consistency of the profiles (|2"0"f is possible only if a takes a special value 
for each coefficient of restitution a. In contrast, the reduced thermal gradient e will remain free 
and so independent of a. 

In terms of the variable s, the temperature profile is 



T(s) = To ( 1 + ) , w = y/2T,,/m. 



(23) 



where To and eo are the temperature and Knudsen number at a reference point s = 0. In 
order to get the dependence of temperature in real space, we need to make use of the precise 
temperature dependence of w. According to the assumption to oc nT 1 / 2 = P T~ 1 / 2 [seeEq. ([XT])], 
Eq. CH]) implies 



yo 

2/o 



— f da' VT(s')/T 
wo Jo 



2v 
3e w 



(T/T ) 



3/2 



where we have called ojq to the value of to at s = 0. Inversion of Eq. (|2"4"|) yields 



T(y) = T 



360^0 

2v 



{y - yo) 



2/3 



Therefore, T 3 / 2 is a linear function of y. Note that e oc dT/dy: 

vo dT 
w T dy ' 



(24) 



(25) 



(26) 



Equations (|2"4"l) - (f2"6")) are only valid if uj oc tiT 1 / 2 (Maxwell particles mimicking hard spheres), 
while Eqs. (|20| -(l23 |) are "universal" in the sense described above. In the remainder of this 
section we will not need to specify the temperature dependence of the collision frequencies. 
Wc now introduce the dimensionless velocity distribution function 



ip(c; e) 



T[s) 
P 



2T{s) 



<I/2 



/(s,v), 



U{8) 



(27) 



Upon writing ip(c; e) we have not included the dependence on the reduced shear rate a since, 
as said before, it is not a free parameter in this solution but it is enslaved to the value of a. 
Since, according to Eq. (12~T)) . all the dependence of / on s occurs through the hydrodynamic 
fields T and u x , we have 

df 8T df , du x Of 



where 



ds ds dT ds du x 
df p / m \( d + 1 )/ 2 dip 



df p / m \ d l 2 
df = T V2T/ 

Taking into account that dc/dT = — 
d 



d + 2 1 9c dip de dip 
2 V+ df'd^ + df~d^ 



ds 



p / m \ (<i+i)/2 
T V2T/ 



\T- X c and de/dT 
d 



\T 1 e, one finally gets 



d 



7T' c + e 7T 
dc de 



, d 
dc x 



p(c;a,e) 



(29) 
(30) 

(31) 
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Consequently, Eq. (TTgl) becomes 

2{d+2) 



_,. )l d_ £ d_ 

dc de) ™ dc x 



(I <>)2f2 d J dCl / da ^ a 1(p ( cf ) ( P( cf i) ~ ^( C Mci)] 
= JWM- (32) 



We consider now the following moments of order k = 2r + £, 

M 3r \t,h{e) = I dcc 2r c^ h c h Mc;e) = (c 2r c^ h c h x ), < h < £. (33) 



By definition, Mo|o,0 = 1 > Mo\i,o = M Q \i.i = 0, and M 2 \q.q = f ■ According to Eq. (1521 . the 
moment equations read 



where 



^2r + £ - 1 - ejjj M 2r \e+i,h+a (2rM 2r - 2 \ i+2 , h+1 + hM 2r \i,h-i) = Jir\e,h, (34) 
J 2r \ Lh = j &cc 2r c l y - h 4j[c\ V , V ] (35) 



are the corresponding collisional moments. As we are considering Maxwell models, J 2r \t,h ar e 
given as bilinear combinations of the form M.2r'W,h'^-2r"\£" ,h" such that 2r' + £' + 2r" + I" = 
2r + I. Therefore, only moments of order equal to or smaller than 2r + I contribute to S 2 r\e,h- 
In particular, from Eqs. (TT21 and ([13")) we have [33] 

Jmo = -c| (36) 

^ = -(l + C)^o, 2 ,, + {| h h Z {) ;\ (37) 
J2\i,h = -(^ + ^qM 2]1 , h , h = 0,l, (38) 



3 /„ ~\ . , 1 



3M 2 |i,o, 


h = 


0. 


M 2 \ 1A , 


h = 


1, 


M 2 \i, , 


h = 


2. 


&M 2 |i,i, 


h = 


3. 



Here, 



5 (i+cja^+jj ,_„• (3<j) 



fs C=l±li^, (40, 
u; a 1 + a 



and we have taken into account that 

^ = i + C, ^ = d ~^4l ^ = §(1 + 0- (41) 

u w a 2 u 2 

Let us check that the hierarchy (f3"4"|) is consistent with solutions of the form 

M 2r \ l<h {e) = 2r jr\? r ^e\ ^ Lh) =0iij+£ = odd, (42) 

3=0 

i.e., the moments M 2 r\i.h( e ) of order 2r + £ > 2 are polynomials in the thermal Knudsen 
number e of degree 2r + £ — 2 and parity £. First note that the moment M 2r \e+i,h m the nrs t 
term on the left-hand side of Eq. (|Mf is a polynomial of degree 2r + £ — 1 and has a parity 
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i 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

2 4 6 8 10 12 14 

&=2r+/ 

Fig. 2. Sketch of the sequence followed in the recursive determination of the coefficients fij' e,h \ All 
the coefficients of the same order k = 2r + £ are represented by a common circle. The zeroth-degree 
coefficients /ig 2r '^'' 1 ' ) correspond to the USF. 



different from that of M.ir\l,h- However, the action of the operator e(2r + £ — 1 — ed/de) restores 
the parity and transforms the degree of the polynomial into 2r + £ — 2. The second term on 
the left-hand side obviously has the same parity and degree as M.2r\i,h- Finally, the condition 
2r' + i! + 2r" + t" = 2r + £ assures that the products M.% r iw ^ M-ir"\l" ,h" are polynomials of 
the same degree and parity as those oi M 2 r\e, li- 
lt is worthwhile noticing that the hierarchy reduces to that of the USF problem for 
IMM when e = [3T], while it reduces to that of the conventional Fourier flow problem for 
elastic Maxwell particles when a = 1 (which, as will be seen below, implies a = 0) [27]. The 
general problem (e ^ 0, a ^ 0) is much more difficult since it combines both momentum and 
energy transport. However, the moment hierarchy can be exactly solved via a recursive scheme. 
The key point is that, while Eq. (|34p includes the moment M.2r\l+i,hi which is a polynomial 

of degree 2r + I — 1, only the coefficients jU j 2r ^ +1 '' 1 J with j < 2r + I — 3 contribute to Eq. 
([34l) . Therefore, one can compute the independent coefficients ^ r ^' h ^ (which are in fact the 

only ones appearing in the USF), then the first-degree coefficients fj^ r '^' h \ and so on. The 
recursive scheme is sketched in Fig. [21 where the arrows indicate the sequence followed in the 
determination of fj,j' e ' h '. The open (closed) circles represent the coefficients associated with 
even (odd) j and £. Notice that the number of coefficients needed to determine a given moment 
■M.2 r \£,h i s finite. In fact, the coefficients represented in Fig. [5] are the ones involved in the 
evaluation of M.2r\l,h f° r & = 2r + £ < 8. In the next Section we follow this scheme to determine 
the pressure tensor (k = 2) and the heat flux (k = 3). 

It must be noted that, although all the moments can in principle be evaluated recursively, the 
explicit form of the reduced distribution function y?(c; e) is not known. However, one could take 
advantage of the knowledge of the first few moments to construct an approximate distribution 
function compatible with them by using, for instance, maximum-entropy arguments. 
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Before closing this section, let us consider the spatial dependence of the dimensional mo- 
ments 



M 2r]i , h (y) = I dvV 2r V«- h V*f( y> v) 



T(y) 



rn 



r+i/2 



M 



2r\i,h 



(43) 



Taking into account Eq. (j4*2"j) and the fact that ^jT{y)e(y) = %/To £ o = const, one can easily get 



M- 



2p 



1 1 \t,h(y) = — e o — 

m \ m 



,,-r \ r+e/2-1 [(2r+i-2)/2] 
Z1 Q \ (2r\t,h) 

/ j L L 2r+e-2-2j 



T(y) 



(44) 



where here [x] denotes the integer part of x and it is understood that 2r + £ > 2. Equation ((31]) 
shows that M 2r ^ y f l (y) is just a polynomial in T(y) of degree [(2r + 1 — 2)/2]. In particular, the 
second- and third-order moments are uniform, the fourth- and fifth-order moments depend on y 
via a linear dependence on T(y), and so on. The spatial uniformity of the third-order moments 
implies that the LTu solution characterized by the hydrodynamic profiles (|20l) corresponds 
indeed to planar Couette flows with uniform heat flux. 



4 Pressure tensor and heat flux 

Let us start considering the moments of second order (pressure tensor) . This will give us the 
dependence of the reduced shear rate a and the rheological properties on dissipation. 

Taking (2r\£,h) = (2|0,0), (0|2,0), (0|2,1), and (0|2,2) in Eq. one gets a coupled set 

of linear equations for o and Mo\2,h {h — 0, 1, 2). The solution is 



\/y(l HI. 



(45) 



M 



0|2,1 



— \ c x c y) — 



M \ 2 ,a = (c 2 y ) 



|2,2 



{cl) 



2(1 + C) 2 

1 



2(1 + 

i + rfC 

2(1 + 0' 



(46) 
(47) 

(48) 



These quantities depend on the coefficient of restitution a through the scaled cooling rate Q 
given by Eq. (|4"0)) . 

Next, we consider the third-order moments. Setting (2r\£,h) = (2|1,0), (2|1, 1), (0|3,0), 
(0 1 3 , 1), (0 1 3 , 2), and (0|3, 3) in Eq. (f3"4"|) . we obtain a set of six linear equations whose solution 
gives the third-order moments in terms of the independent terms of the fourth-order moments 
and Here we only display the explicit forms for the two third-order moments defining the x 
and y components of the heat flux. They are 



M 



2|1,0 



(0 



(c 2 c„ 



2de 



M{2d -2 + 3dC)CMo°' 4 ' 0) + 8dCAio°' 4 ' 2) + (18d - 18 + 19<*OMo 



(2|2,0) 



6\/2dC (2d-2 + 3dOMo 



(0|4,1) 



(212,1) 

Mo 



(49) 
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M 2]1<1 (e) = (c 2 c x ) = |||^2dC[M7d - 2 + 9dC)C^° |4 ' 0) + 6(6d - 6 + 5d()^ 4 > 2) 
+3(9d + 6 + 17d()fi { o lm ] - 12d({7d -2 + 9dCV 0|44) 
-9(6d-6 + 5dC)^ 2|24) \, (50) 



where 

X = 36(d-l) 2 -d(j6-56d-9d?)C (51) 

Thus, in order to complete the determination of .M2|i,o and A^2|i,i we need to know the a- 

dependence of the coefficients /z °' 4 ' G \ /z °' 4 ' 2 ', /i 2 ' 2,0 ', /i ' 4,1 ', and ^ 2,1 \ These coefficients 
can be obtained from Eq. (fMf by taking e = 0. This was in fact already done in Ref. [3T] in 
the context of the solution of the USF for IMM. The final expressions for the three-dimensional 
case (d — 3) are 

■M 2 |i,o(e) 135(1 + a) 2 A(a) 



8(4 - a) 2 (829 - 162a - 91a 2 ) C(a) ' 



■M 2 |i,i(e) _ 27V5/2(1 + a)Vl - a 2 B{a) 
e ~ 4(4 - a) 2 (829 - 162a - 91a 2 ) C(a) ' 



(52) 
(53) 



21. 



where 

26 

A(a) =5^(1 -a) 1 , S(a) = ^B i (l-a) i , C(a) = ^ C,(l - a)'. (54) 

i=0 i=0 i=0 

The numerical values of the coefficients Ai , £?i , and Ci are given in Table [TJ 

Once the pressure tensor and the heat flux have been determined, it is convenient to define 
generalized transport coefficients. The relevant elements of the pressure tensor are characterized 
by an effective (reduced) shear viscosity 77* and two independent (reduced) normal stresses 9 X 
and 8 y defined as 

P xy = -V* — ^, (55) 

^ = e x , ^ = e y . (56) 
p p 

If d > 3, the other d — 2 normal stress coefficients 9 Z = ■ ■ ■ = 9d are related to 9 X and 9 y 
by 9 Z = (d — 6 X — 9 y )/(d — 2). The non-zero components of the heat flux q y and q x define an 
effective (reduced) thermal conductivity A* and a (reduced) cross coefficient c/>* , respectively, 

by 

., d(d + 2) p dT _^ d(d + 2) p dT 
% ~ A 2(rf-l)m^s^' qx ~ 2(d-l)m^ s dy- 

From Eqs. (g5]h(|IH]) we easily get 



2 



1 du x f 1 + a\ 



d+l-a d+2 

= ^2d— V^ (1 - a) ' (58) 

^^(^) 2 ' (59) 

W = L+l-a ' U*) = d-{d-l)e y {a), (60) 



Will be inserted by the editor 
Table 1. Numerical values of the coefficients Ai, Bi, and d defined in Eq. i|54fl . 
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i Ai Bi Ci 






341218675233128448 


2985663408289873920 


2369574133563392 


1 


2482013713849384960 


20450944454204522496 


7160814742552576 


2 


7969339170892868608 


58568421476117803264 


5980987758379200 


3 


14758809962784485248 


96685932126035252672 


-1622483991362528 


4 


17947525965827798144 


107159766670778436160 


-6127404964857736 


5 


15792893280853221008 


88222721223517916056 


-5323957544763312 


6 


10848492042169821200 


57714899464650979880 


-3091956803616730 


7 


6083138827020141820 


31096845798718075642 


-1463253462204164 


8 


2852642631673202332 


14152867760388474114 


-595099812423405 


9 


1160168028460008358 


5649360969446546811 


-211065154772581 


10 


424875366505855228 


2030043005756227200 


-63304902022304 


11 


137908868657776387 


642139536452566558 


-15813143468850 


12 


37996059299305384 


172739779638915362 


-3363854349023 


13 


8835881256127730 


39251245670664391 


-521942377705 


14 


1676141413944428 


7120337713266190 


-23351656362 


1 r. 


iy41zolooUoOzl ( 


( z44 f oU4oloz4oo 


1 A T A A Q QQQflO 

14 ( 44ooyoUz 


16 


-14258415439016 


-96248350169696 


7223865310 


17 


-15161725120946 


-74193692948568 


2126089824 


18 


-5457337268162 


-25236538926270 


429672915 


19 


-1406011921458 


-6243932259627 


78316371 


20 


-274206254448 


-1200674292216 


11377737 


21 


-47059936911 


-201709797345 


1142559 


22 


-6553863486 


-27252518922 


141858 


23 


-661870251 


-2792296044 


5508 


24 


-75802878 


-300657636 


648 


25 


-3147012 


-13199760 
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-329832 


-1263600 






where use has been made of Eq. (HOI) . In Eq. (I55|) . a is the shear rate scaled with the collision 
frequency z^hs- The second equality of Eq. (joT))) implies that 6 Z = 9 y (if d > 3). Note that, from 
Eqs. CE7j), flM]), (EH), and dSHJ), we have 



V (a)a (a) = - , (61) 

2 i^HS 

in consistency with d y q y =0 in Eq. ([1}. 
The heat flux coefficients are given by 

16(rf-l) -M 2lhQ 16(d-l) -M 2 \ U1 

X{a)= d(d + 2)(l + ar e ' * {a) = d(d + 2)(l + ay e ' (62) 



More explicitly, for d = 3, 



[ ' (4-a) 2 (829 - 162a-91a 2 )C(a)' 1 ' 



,h*( ) = 7275/271^ Bja) 

9 1 ' 5(1 + a) (4 - a) 2 (829 - 162a - 91a 2 ) C(a) ' 1 J 

where the functions A(a), B(a), and C(a) are given by Eq. (l54l) and Tablc[TJ 



5 Comparison with computer simulations 



The exact solution to the BE @ and ([3]) obtained by the moment method in Sec. |3] defines a 
normal or hydrodynamic solution where the spatial dependence of all the moments only occurs 
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0.5 0.6 0.7 0.8 0.9 1.0 

a 

Fig. 3. Plot of the reduced shear rate a (top panel), the reduced shear viscosity rf (middle panel), and 
the reduced normal stresses 6 X (bottom panel, upper curve) and 6 y (bottom panel, lower curve) for a 
three-dimensional system. The lines are the theoretical results (|58p - (|60p . while the symbols correspond 
to DSMC results using two different applied thermal gradients: AT/T w d — 2 (x) and AT/T wc i = 5 
(+)• 



through the hydrodynamic fields n(y), u x (y), and T(y), as shown by Eqs. ([27|). and (flU) . 
In fact, this solution is free from boundary-layer effects. 

The important point is whether or not the exact solution derived here actually describes 
the steady state reached by the system, in the bulk domain, when subject to realistic boundary 
conditions and for arbitrary initial conditions. To confirm this expectation, we have obtained 
the numerical solution of Eqs. @ and ^ by means of the DSMC method [TU] for a three- 
dimensional id = 3) system. 

The DSMC algorithm consists of two basic steps, both iterated each time interval St: (a) a 
free streaming step, in which all particles are drifted by vSt (v being the particle velocity), and 
(b) a collision step, in which vbt collision pairs are selected among neighbor particles, where v 
is the characteristic collision frequency appearing in Eq. and given by Eq. (1111) to mimic 
hard spheres. Additionally, for a bounded system like ours (see Fig. [1]), the particle velocities 
are updated during step (a) if they eventually touch the boundaries, which in our case are 
thermal walls moving with constant and opposite velocities (see more details on the description 
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^ I i i i | i i i i | i i i i | i i i i | i i i i | i i i i | i i i i f-^ 

s 

S 

4: X x : 



* x 




0.70 0.80 0.90 1.00 

a 

Fig. 4. Plot of the reduced heat flux coefficients A* (solid line and filled triangles) and 0* (dashed 
line and open triangles) for a three-dimensional system. The lines are the theoretical results (|63p and 
(|64p. while the symbols correspond to DSMC results using two different applied thermal gradients 
(AT/T wc i = 2 and 5). 



of these boundary conditions elsewhere [3D] )■ The difference between DSMC for the Maxwell 
model and DSMC for hard spheres lies in the collision step. Since now the collision probability 
is independent of the relative velocities of the colliding particles [see Eq. ([3])], the collision pairs 
are randomly selected with equiprobability. Once the particles are selected for collision their 
velocities are updated following the same collision rule as for hard spheres [see Eq. ([5])]. The 
number of particles N and the time step St need to be sufficiently large and small, respectively 
[ID] . We use here N = 2 x 10 5 and St = 3 x 10 _3 z/ _1 , where v is given by Eq. (fTT]) with n — > h 
and T — > T w d, n and T w d being the average density and the cold wall temperature, respectively. 
In addition, the system needs to be split into small enough cells [TDJ and we chose layers of 
width Sy = 2 x 10 _2 A m f p , A m f p = (-K^f^fuj 2 )^ 1 being the characteristic mean free path. 

In order to find the states with uniform heat flux in the bulk, we run series of simulations with 
varying wall relative speeds U wu — U w d for fixed wall temperature difference AT — T wu — T w d 
and distance h. This results in series of varying shear rate |8l9j . Since the theoretical values 
of the reduced shear rate a(a) for IMM are very close to those for IHS, we have used values 
of h and AT in the same ranges as those used in Ref. [DJ: h/X m fp ~ 10-20, AT/T w d = 1-20; 
the appropriate range of U wu — U w d depends on the value of a [9]. In order to obtain steady 
states we need to wait for a sufficiently long computing time. In the case of Maxwell particles 
this time turns out to be much longer than for hard spheres since the third-order moments 
are coupled to the fourth-order moments of the USF, which have very long relaxation times 
[21j . Also, for better averaging the steady values of the hydrodynamic profiles and transport 
coefficients, we average twice: in time, by means of repeated measurements in the steady state 
at different uncorrelated instants; and in space, by averaging the small simulation cells to larger 
hydrodynamic cells, analogously to how we proceeded in a former work [30] . The rest of technical 
details are the same as in [3D] , 

Figure [3J shows the reduced shear rate a, the reduced shear viscosity 77* , and the reduced 
normal stress coefficients X and 9 y as functions of the coefficient of restitution a for a three- 
dimensional system. The DSMC results obtained by using two different applied thermal gradi- 
ents are also plotted. Computer simulations show that, for each value of a, the values of a{pt), 
77* (a) and 9i(a) are insensitive to the choice of the thermal gradient. This confirms the inde- 
pendence of the rheological properties and the shear rate on the thermal gradient, including 
the case of no thermal gradient (USF). We also observe that the agreement between theory and 
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simulation is excellent, even for quite strong values of dissipation. Regarding the influence of 
dissipation on transport properties, we see that the shear viscosity is significantly affected by 
dissipation since the deviation of rj* from its elastic value (rj* — 1) is quite relevant even for 
moderate dissipation (say a w 0.9). It is also important to notice that such a deviation from 
the value in the elastic limit has a sign opposite to that of the NS shear viscosity [3D]. With 
respect to the normal stress coefficients, Fig. [3J shows that the anisotropy induced by the shear 
flow is more important in the direction parallel to the flow velocity (x direction) than in the 
direction parallel to the thermal gradient (y direction). 

The a-dependence of the (reduced) heat flux coefficients A* and </>* is plotted in Fig. @] 
As said before, given that the transient regime needed to reach steady-state values of these 
coefficients is longer than that of the rheological properties (especially for strong dissipation), 
the simulation data reported here (a > 0.7) cover a range of values of a smaller than the one 
plotted in Fig. [3J However, a > 0.7 is still of course a relevant range of values of dissipation 
since most of the experiments carried out for granular fluids are included in this region. Fig- 
ure 2] clearly shows again that the heat transport is clearly affected by collisional dissipation, 
especially in the case of the cross-coefficient <jf . As happens for the rheological transport co- 
efficients, the theoretical results compare very well with simulations in all the range of values 
of a studied. Interestingly, it can be be observed that the streamwise component q x becomes 
larger in magnitude than the crosswise component q y for a < 0.9, what represents a strong 
far-from-equilibrium effect. 



6 Summary and Discussion 

Exact solutions in kinetic theory for far-from-equilibrium states are exceedingly rare. The dif- 
ficulties increase in the case of a granular gas of smooth IHS since a new parameter (the 
coefficient of normal restitution a) is introduced to account for the dissipative character of col- 
lisions. Usually, in order to get explicit results from the BE one considers approximate methods 
such as the leading order in a Sonine polynomial expansion of the velocity distribution func- 
tion or the well-known Grad's 13-moment method. Although in many of the cases the results 
obtained from the above analytical methods compare reasonably well with numerical solutions 
of the BE, the lack of exact solutions motivates the search for alternative routes. One possibil- 
ity consists of replacing the detailed Boltzmann collision operator by a simpler collision term 
(e.g., the Bhatnagar-Gross-Krook model kinetic equation), that otherwise retains the most 
relevant features of the true collision operator. Another different approach consists of retaining 
the mathematical structure of the Boltzmann operator, but considering a simplified collision 
model, such as the IMM. As in the case of elastic collisions [26|27| . the collision rate of this 
model is assumed to be independent of the relative velocity of the two particles that are about 
to collide. 

The simplicity of IMM allows us, for instance, to get exactly the collisional moments of the 
BE without the explicit knowledge of the velocity distribution function [23] . Their knowledge 
has opened up the possibility of exactly determining the nonlinear transport coefficients in 
some far-from-equilibrium situations, such as the USF problem |19l21j . In this paper, an exact 
solution for the class of steady Couette flows defined by a uniform heat flux has been found. 
This solution is characterized by a constant pressure p = nT and by linear profiles of the as- 
component of the flow velocity u x and the temperature T with respect to a conveniently scaled 
space variable. As a consequence, a linear relationship between T and u x holds, so that we 
refer to this class of steady Couette flows as the LTu class. Moreover, as in the USF state, the 
condition of stationarity imposes a relationship between the shear rate and the coefficient of 
restitution, so that the (reduced) thermal gradient e is the only free nonequilibrium parameter 
of the problem. In particular, this steady Couette flow reduces to the USF when e = (with 
a^l), while it reduces to the conventional Fourier flow problem when a = 1 (with e ^ 0). A 
previous study of the LTu flow problem has been recently carried out by the present authors 
for IHS [9 by considering three complementary routes: an approximate solution from Grad's 
13-moment method, DSMC simulations, and molecular dynamics simulations. 
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A careful inspection of the moment hierarchy (|34j) (seen as an infinite set of first-order 
differential equations with respect to e) shows that it can be exactly solved via a recursive 
scheme. The solution is characterized by reduced moments of the form (|42p . i.e., moments of 
order k are polynomials of degree k — 2 in the (reduced) thermal gradient e. The coefficients 
of these polynomials are nonlinear functions of the coefficient of restitution a that can be 
recursively solved starting from the USF values (where the hierarchical character of the moment 
equations is broken down) and following the scheme sketched in Fig. [2] In particular, the 
elements of the pressure tensor are independent of e and their expressions coincide with 
those of the USF [3T] . The above elements define an effective (reduced) shear viscosity rf and 
two independent (reduced) normal stresses X and 9 y . The heat flux q is linear in e and has 
two nonzero components: one component q y parallel to the thermal gradient and a component 
q x normal to the thermal gradient. In reduced units, the first component defines an effective 
thermal conductivity A* while the latter component defines a cross-coefficient 0* that is absent 
in the NS description. The evaluation of the explicit dependence of A* and <fi* on a has been one 
of the main goals of this paper. This dependence has been obtained thanks to the knowledge 
of the fourth-degree velocity moments of IMM in the USF [3T] . As expected, the results show 
that the influence of collisional dissipation on the transport coefficients -q* , 0j, A*, and <f>* is 
quite significant (see Figs. [3] and E]). This is especially apparent in the case of cf>* (which is zero 
in the elastic limit) since its magnitude becomes even larger than that of the effective thermal 
conductivity A* for moderate values of dissipation. This is a signal of the strong non-Newtonian 
behavior of the LTu flows for the granular gas. 

In order to validate the theoretical results derived from the moment method, computer 
simulations based on the DSMC method have been carried out for the same interaction model 
in the three-dimensional case. Given that the exact solution derived in Section [4] is free from 
boundary-layer effects, the important point to elucidate via comparison with computer simu- 
lations is whether or not this analytical solution actually describes the steady state reached 
by the gas in the bulk domain, when the gas is subject to realistic boundary conditions and 
for arbitrary initial conditions. The good agreement between theory and simulation found here 
confirms the hydrodynamic profiles and the nonlinear transport properties predicted by the 
analytical solution, even for quite strong values of dissipation. 

This work has been supported by the Ministerio de Education y Ciencia (Spain) through grant No. 
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